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CN ; Abstract 

The density of states method is applied for lattice QCD at a finite isospin density. The advantage 
1 of this method is that one can easily obtain results for various values of parameters ( quark mass, 

coupling constant and the number of flavors ). We compare results for the chiral condensate and the 
quark number density with those from the R- algorithm and find that they are in good agreement. 
By calculating the chiral condensate we obtain information on the phase structure for various 
■ quark flavors and isospin chemical potentials. We also show results for the chiral condensate at two 

^ \ different quark masses and at two different isospin densities which are not easily obtainable in the 

conventional Monte Carlo method. 
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1 Introduction 

CO . ... 

The lattice QCD Monte Carlo technique has been a useful tool for clarifying the non-perturbative 
aspects of QCD at the zero baryon chemical potential ( \ib = ) for both finite and zero temperatures. 
At nonzero hb, however, due to the sign problem, the standard importance sampling method fails. 
The Glasgow group 1. attempted to extract information on hb at low temperature by a method 
based on reweighting but this method does not work at low temperatures due to the overlap problem. 
Recently it has been realized that at low density and finite temperature the sign problem may not 
be a serious numerical difficulty and a variety of approaches ( multi-reweighting, Taylor expansion, 
imaginary chemical potential) 03 IU 03 E] have been applied to study QCD at low density and 
finite temperature (See 03 03 EH for recent reviews). For larger \ib we are still lacking an efficient 
numerical method. 

In contrast to the finite \ib case, QCD at isospin density has no sign problem and can be simulated 
by lattice Monte Carlo methods 11 . There is an expectation that at small isospin density the system 
may resemble that at small baryon density [B| 112] . Recent studies have shown that the phase diagram 
at small isospin density is similar to that at small jhr |131IT3| . On the other hand, a difference between 
the isospin and baryon density is seen in susceptibilities of number density and derivatives of meson 
mass P ITS]. 

In lattice Monte Carlo simulations, usually importance sampling is used. This approach have been 
proven to be very efficient by many studies. Another approach one may take is the density of states 
(DOS) method. The difference between importance sampling and the DOS method is as follows. 
In importance sampling a simulation is performed at a fixed parameter set. To explore different 
parameters one must perform independent simulations 1 . On the other hand, in the DOS method one 
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J If the parameter region to be explored is narrow, then one may apply the reweighting method[T7)j which does not 
need independent simulations. 
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first determines the DOS of suitable observables and then results are obtained by performing one or 
a few dimensional integral. To illustrate how the DOS method works, let us consider a gauge model. 

The partition function for this model is given by Z = J [dU]exp(—/3S g [U]), where S g [U] is a gauge 
action and (3 is a coupling constant. If we define the DOS by n(E) = J [dU]5(E — S g [U]), then the 

partition function is rewritten as Z = J dEn(E) exp(— j3E) . If we consider the average value of the 

gauge action, it is given by < E >= J dEn(E)E exp(—PE)/Z . This is a one-dimensional integral of 
E and once we obtain n{E), we can evaluate this easily at any (3. No extra simulation is needed to 
calculate it at various /3, which is considered to be an advantage of this method. 

The DOS method has been applied for gauge theories: Z(n), SU(2) and SU(3) models |T7j. Although 
the inclusion of dynamical fermions is computationally difficult, in Ref. the DOS method was 
applied for QED, where the method is called the microcanonical fermionic average method. Recently, 
LuoilQ] extended the idea of the DOS method to QCD and emphasized that once the eigenvalues 
of the Dirac operator and the DOS of the plaquette energy are determined, one can evaluate the 
thermodynamic quantities derived from the partition function at any quark mass and flavor. 

In this study we apply the DOS method for QCD at isospin density on a 4 4 lattice and demonstrate 
how the DOS method works. Since the isospin system has no sign problem, results from the DOS 
method can be compared with those from the standard method. In Sec. 2 we give general formulas 
for forming the DOS including dynamical fermions. In Sec. 3 we give simulation details. Results are 
presented in Sec. 4. We summarize our results in Sec. 5. 



2 Density of States Method 
2.1 General Formulas 

In lattice QCD Monte Carlo simulations, usually we aim at obtaining the average values of observables: 

(O) = |y [d[/]0[C/]detA(m (? ,/x) Af // 4 exp(-/?^[[/]), (1) 

where Nf is the number of flavors and A(m q , //) is assumed to be a staggered fermion matrix at quark 
mass m q and at chemical potential /i: 

A(m q , fi)i t j = m q 6 it j + - Vi,v( u i,v$i,j-v ~ Ui_„ jU 6ij +u ) (2) 

!/=l,2,3 

S g [U] is the standard Wilson gauge action, 

S 9 [U] = -^J2 ReTr ( U p)> (3) 

Jv c p 

where U p is the plaquette and N c = 3 for SU(3) gauge theory. For Wilson fermions, Nf/4 in eq.Q 
should be replaced with Nf. The standard means of dealing with this multi-dimensional integral is to 
do the importance sampling, i.e. configurations are generated with a measure ~ dll det A N ^ A e~^ Sg 
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and (O) is given by an average over the configurations as (O) ~ J2i=i 0[Ui], where N is the number 
of configurations. 

The DOS method rewrites the partition function and reduces it to a few dimensional integral. Let 
us define the DOS with parameters Ei (i = 1, . . . , k) as 

n{E 1: E 2 ,..., E k ) = J [dU]g{U)IkS(Ei - Xi (U)), (4) 

where Xi(U) is an operator associated with E{ and g(U) is introduced to generalize the DOS further 
with g(U) ^ 1. Using the DOS, we obtain 

(°) = I [^idE i ]n(E 1 ,E 2 ,...,E k ){OdetA N f/ 4 e- l3 ^/g(U)) EuEl _ Ek , (5) 

where 

Z n = J [U i dE i ]n(E l ,E 2 , . . .,E k )(detA N f/ 4 e-^/g(U)) EuE2 ,...,E k . (6) 

Here, (A) EltE2 ,...,E k stands for an average of A on configurations generated with the measure [dU]g(U)TLi5(Ei — 
Xi(U)), or one can write it as 

(A) El , E2 ,..., Ek = * — / dUg(U)Ui6^Ei - Xl {U))A(U). (7) 

n(Ei,E 2 , ... ,E k ) J 

For g(U) = 1, (A)e 1! e 2 —,E k becomes a microcanonical average at fixed E\,E 2 , . . . ,E k . 



2.2 g(U) = l 

Here we give the formulas used in the present study. We consider the case of one parameter and choose 
x(U) = S g [U\. Although there are many possibilities to choose x(U), this choice with g(U) = 1 is the 
one used in Ref.fH?] and is useful for our purpose. Setting g(U) = 1, we obtain 



n(E) = J[dU]5(6VE-S g [U]), (8) 

(0) = ^- f dEn(E)e- p6VE (OdetA(fi) N f/ 4 ) E , (9) 

and 

Z n = J dEn(E)e-^ VE (detA(p) N f/ 4 ) E , (10) 

where V is the number of lattice sites and E is the plaquette energy. These are the basic formulas 
used in our study. Our task includes three numerical calculations: (i)n(E), (ii) (det A(fj,) Nf ^' i ) E and 
(iii) (OdetA(fi) N ^ 4 ) E . Luo argued that if one stores the eigenvalues of the fermion matrix for 
all configurations, then one can evaluate (det A(n) Nf l 4 ) E for any quark mass and flavor. Let Aj(/x) 
be the i-th eigenvalue of the massless fermion matrix A(m q = 0) on a configuration with E. Then we 
obtain 

(det A{ii) N f' 4 ) E = (( ft (W + m q )) N f/ 4 ) E . (11) 

i 

Since (iii) contains O, in general it is not calculable for any quark mass and flavor. However the chiral 
condensate {tpip) which is obtained as the trace of the inverse fermion matrix can also be given with 
the eigenvalues. Namely, for ("/pip) we obtain 

1 1 1 NcV 

((-TrA(^)- 1 ) det A(fi) N ^ 4 ) E = (- ]T ^ + ( II + m q )) N f' 4 ) E . (12) 
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2.3 g(U)^l 



Depending on the purpose, one may take g(U) ^ 1. In Ref.[2T)] the DOS for the complex phase of the 
fermion determinant was calculated. The DOS was defined as 

n(E) = J [dU]g(U)S(E - 6(U)), (13) 

where g{U) is set to | det /\\e~^ Sa<yU \ and 9(U) is the complex phase of the determinant. Then the 
expectation value of O is given by 

(O) = r dE(0) E n{E)e iE / f dEn(E)e iE . (14) 

J — TV J —7V 

In Ref.[2~T] the DOS for the number density was used and g(U) is also set to | det A|e _ ^ 59 ^^. 

In these definitions parameters (3 and m q are absorbed in n(E) and we cannot vary (3 and m q . 
Therefore, in this study we do not use these definitions with g(U) ^ 1. 



3 Simulations 

For the case of Nf = 2, we have two chemical potentials fj, u and /i^ associated with the u and d quarks 
respectively. If we introduce the isospin chemical potential uj as [ii = [i u = — we obtain 

[detA^A^)] 1 / 4 = [det A( W )A(- W )] 1 / 4 = | det A^)] 1 / 2 , (15) 

which is positive in the Monte Carlo measure. Therefore, for this \xi the sign problem is absent and 
one can perform Monte Carlo simulations. The natural generalization of the isospin chemical potential 
to N f ^ 2 is that one takes det A^/)^/ 4 = | det A^/)^// 4 . 

We follow the implementation developed in Ref.49 a . The DOS n(E) in eq.(jHJ) can be obtained 
using the quenched data as 

- = 6 f E d E'f3(E') + const. (16) 

V Jo 

First we make quenched simulations on a 4 4 lattice and determine the coupling constant /3(E) as a 
function of the plaquette E ( Figure^). Then we integrate (3(E), numerically interpolating the data 
according to the trapezoidal rule. Figure |5] shows the result of — lnn(E)/V. 

The most time consuming part of our method is the calculations of the microcanonical averages 
(det A(m) N f^) E and (O det A(/i / ) 7V // 4 } i3 function of E, which contain the eigenvalue calculations. 
In order to generate configurations at E we use the over-relaxation method [22} I19j. Starting from a 
configuration with E, we have made the over-relaxation update and saved 100 configurations at each 
E. Each configuration is separated by 100 over-relaxation updates. About 30 values of E are chosen 
in E E [0.05, 0.95]. For each configuration we calculate the eigenvalues of the fermion matrix and store 
them. Using these eigenvalues we can evaluate (det A(m) N f^) E as in eq. (fTTj) for any quark mass and 
flavor. Examples of (det A (fj,j) Nf 7/4 } E are shown in Figure El Eq. (|12j) is evaluated similarly. These 
data are interpolated by polynomials when we perform eas.(|8])-([10|). 
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Quenched simulations on a 4 lattice 




Figure 1: Plaquette energy as a function of (3. 




<E> 

Figure 2: — In n(E)/V as a function of plaquette energy E. 

4 Results 

Figures 0] and |S] compare results of {^ip) between the DOS method and the R- algorithm |23| . We see a 
good agreement between them in a wide range of /3. A small difference is seen in the phase transition 
region where the {ipip) changes rapidly. In such a region we probably need careful analysis. 

Figure El shows (ipifj) for different Nf at m q = 0.025 and [ii = 0.2(left), and at m q = 0.05 and 
[ii = 0.25(right) as functions of [3. One can see that the critical coupling where the phase transition 
occurs decreases as Nf increases. Figure [7| shows (t/jip) for various /i/. In the low-temperature phase 
{4>ip) decreases as m increases. On the other hand, in the high-temperature phase no visible difference 
can be seen among (ipip). Since the DOS method can explore various parameter space easily it is 
considered to be useful for exploring a wide parameter space and for seeing a rough phase diagram. 

In the DOS method we can take various combinations of parameters. Let us consider the case of 
Nf = 1 + 1 with non-degenerate quark masses mi and mi- In this case we must calculate the following 
microcanonical averages: 

(IdetAMI^IdetAMI^/ 4 )^ (17) 
WK=i )2 )| det A( mi )\ N f^\ det A(m 2 )\ N f/ 4 ) E . (18) 
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Figure 3: Micro canonical average of det A(m) N f' A at m = 0.2 and m q = 0.025 as a function of E. 
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Figure 4: (i/>^>) for Nf = 2 at m q = 0.05 and at /i/ = O.O(left) and 0.2(right). Circles are from the 
DOS method. The results from the R-algorithm are shown with squares. 



Since the eigenvalues are stored, it is easy to calculate these microcanonical averages. On the other 
hand, in the conventional algorithm such as R-algorithm one needs a differently implemented program 
to simulate Nf = 1 + 1. Such a program can be implemented but may become intricate. 

Figure |H1 shows {ipi^) for Nf = 1 + 1 with different quark masses ( mi = 0.05 and m.2 = 0.025 ) 
and at jxi = 0.2. Since the two quarks have different masses, we plot two {tpip(m q )) for mi and mi. 
Similarly, we can also consider non-degenerate isospin chemical potentials ( \x\ and M2 )■ In this case 
we calculate 

(| det A( Ml )|"// 4 | det A(n2)\ N f/%, (19) 

(^(w=i, 2 )|detA( / xi)| 7V // 4 |detA(/i 2 )| Ar // 4 ) s . (20) 

Figure shows (ifjip) for Nf = 1 + 1 with different chemical potentials ( \i\ = 0.2 and \ii = 0.3 ) at 
m q = 0.025. 

Next we show the results of the number density defined by 
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Figure 5: (ijjip) for N f = 4 at m = 0.2 and at m g = 0.05(left) and 0.025(right). 
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Figure 6: (V^) for different TV/, (left) m q = 0.025 and \ii = 0.2. (right) m 9 = 0.05 and m = 0.25. 



Since this can not be evaluated with the eigenvalues of A(/i) only, in order to obtain (0\ det A(fj,)\ N f^ 4: ) e 
for the number density we must calculate the number density on each configuration. 

Usually Eq. (|21|) is evaluated by the noise method. Using noise vectors Ri having the property 
(RjRj) = 5ij, the trace calculation can be replaced by 



Tr 



1 dA(fi) 



ifnf 1 ^fr) 

N<H 1 A{ij,) dfi 



Ri, 



(22) 



i=i 



where is the number of noise vectors. Although this is true for — > oo, we find that the convergence 
of the noise method is extremely slow on a 4 4 lattice. A similar result was reported in Ref . |241 . Figure 
HOI shows the convergence measured as relative error as a function of the number of Z2 noise vectors [25], 
Typically, 0(1000) noise vectors are needed to have a reasonable value. In this analysis, instead of using 
the noise method we calculated the number density exactly by calculating each diagonal element of 

- - — . In general such calculations are computationally costly and should be avoided. However, 
A(/x) dfi 

our lattice size is sufficiently small to perform the exact calculation. Thus, here we have adopted the 
exact calculation. 
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Figure 8: (ipip(m q )) for Nf = 1 + 1 at /// = 0.2 with different quark masses, m q = 0.05 and 0.025. 



Figure ITT1 shows the number density as a function of fij. The results from the R-algorithm are also 
plotted in the figure. The two results are in good agreement. 

5 Summary and outlook 

We have given the general formulas of the DOS method including dynamical fermions. The case of 
"g(U) = 1" corresponds to that used by Luojl9j. Based on the implementation by Luo, we have 
calculated {ipip) and the number density at finite isospin densities by the DOS method and made a 
comparison with results from the R-algorithm. The two results were found to be in good agreement. 
The DOS method can be applied for various combinations of parameters. We have calculated (tpip) for 
various values of Nf and /x/, and also for non-degenerate quark masses and for different isospin chemical 
potentials. Especially it is emphasized that for non-degenerate quark masses and for different isospin 
chemical potentials it is not easy to perform Monte Carlo simulations by the conventional algorithm 
as the R-algorithm but these calculations are easily performed in the DOS method by keeping all 
eigenvalues of the fermion matrix. 

The limitation of the DOS method may appear on a large lattice. The measurement of the DOS 
method is done in the microcanonical ensemble and (det A) is treated as an observable. Since the 
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Figure 9: (ipip(fii)) for iVy = 1 + 1 at m g = 0.025 with different \i\ = 0.2 and 0.3 as a function of (5. 



microcanonical and full ensembles are very different, (det A) is expected to fluctuate largely as the 
volume of the system increases, which may limit the available lattice to a small one. 

In principle, the DOS method can be applied for QCD with a baryon chemical potential. How- 
ever, there still remains the sign problem. For example, for Nf = 4 we have (det A(^lb))e = 
(| det A(nB)\e ie )E and if e l6 fluctuates significantly, which is expected to occur for larger /^b, one cannot 
obtain meaningful values for the microcanonical average. We attempted to calculate (| det A(/j,B)\e ie ) e 
but for fj-B > 0.2 we could not obtain statistically meaningful values 2 . On the other hand, for [Ib < 0.2 
the results were stable but there is no visible difference between the results of (ipip) at \xi and at hb', 
for [ib < 0.2 the effect of the complex phase is small and thus we do not see any difference. This is 
consistent with the results of Refs. |2lfl 12*7] . 

In this study for each value of \i we calculated the eigenvalues. However, if we use a matrix reduced 
to two time slices 28, E|, we can calculate the determinant at any [i. The determinant for this case is 
expressed as 

det A(n) = det(P + e Ntfl ) x e 3 ^, (23) 

where P is a /i-independent matrix and Nt is the number of lattice sites in the time direction. Keeping 
the eigenvalues of P, we can calculate the determinant at any \i. In this case, however, the tradeoff is 
that the determinant is not calculable at any m q . Thus, if one wishes to obtain the behavior with \i 
at fixed m q this reduction method may be useful. 
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